!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! On the Optimal Design of Transfers and Income-Tax Progressivity
! Grid search for m-theta-line: No Pareto with no phase-out
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program Main_PreTesting

    use omp_lib
    use Globals
    use Toolbox
    use ModuleSolveEqm
    use ModuleQTRANSITION

    implicit none 
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Local variable declaration
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    ! For submission of many programs
    character(len=5)    :: seq_no_s     ! read input from Linux directly
    integer             :: seq_no       ! sequence number of program 
    integer             :: no_assign    ! number of steady states assigned to each program
    integer             :: loop_start   ! index of first steady state computed by this program
    integer             :: loop_end     ! index of last steady state computed by this program
    logical             :: exist        ! to check whether file already exists

    ! For optimization
    integer, parameter :: N = 2                 ! number of parameters to be optimized
    integer, parameter  :: nmt = 21             ! number of points for transfer level
    integer, parameter :: ngamma = 7            ! number of points for tax progressivity
    integer, parameter :: II = ngamma*nmt       ! number of grid points
    integer, parameter  :: no_prog = ngamma     ! number of programs to be run (II needs to be multiple of no_prog)
    integer :: ia, ian                          ! counter
    real(8) :: paramSOB(II, N)                  ! parameter matrix
    real(8) :: fvalSOB(II)                      ! welfare at the equilibria 
    real(8) :: fvalSOB_ss(II)                   ! welfare only for steady states
    integer, parameter :: N_info = 20           ! number of eqm results to be stored
    integer ::  isave                           ! counter for saving
    real(8) :: info(II,N_info)                  ! matrix of results over grid
    real(8) :: mt_grid(nmt)                     ! grid for transfer level
    real(8) :: gamma_grid(ngamma)               ! grid for tax progressivity
    integer :: ibn, icn, idn, ix                ! counter
    real(8) :: welfare                          ! welfare associated with equilibrium
    real(8) :: param_vec_init(3), param_vec_lr(3)   ! vector of tax function parameters: input for function solving for equilibrium given taxes
    
    ! Errors 
    real(8) ::  errors_vec_Sob(II,4)             ! errors for each Sobol point
    real(8) ::  errors_vec_trans(II,2,T_TR)      ! errors for transition
    real(8) :: error_scalar_trans(II,2)          ! maxval of errors_vec_trans in transition
    
    ! Calibrated parameters
    real(8) :: ParamCalib(5), eqm_results(5)    ! calibrated parameters
    
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Read sequence number of this program
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    if ( COMMAND_ARGUMENT_COUNT() /= 1) then                           
        print *, 'ERROR: A command line argument specifying the parameter is required'
        STOP
    endif 
    print *, 'about to read'
    call get_command_argument(1, seq_no_s)
    print *, 'done reading'
    read(seq_no_s,*) seq_no              
    print *, 'done writing seq_no'            
    print *, 'seq_no    = ', seq_no
    print *, 'seq_no_s   = ', trim(seq_no_s)
    
    
    ! Initialize info file storing results
    if (seq_no == 1) then
                
        info=0
        open(1,  file = 'Output/info.txt',      status = 'unknown')
        write(1, '(*(f18.8))') ( Info(1,isave), isave = 1, N_info )
        close(1)
    
    endif

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Assign calibrated parameters
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    ! Choose normalization version
    norm_version = 1            ! indicator for choice of normalization: 1) mean income; 2) median income; 3) mean income of the third quintile
    
    ! Calibrated parameters
    open(11,  file = 'Input/ParamCalib.txt',      status = 'unknown')
    read(11, *) ParamCalib
    close(11)
    B           = ParamCalib(2)             ! labor disutility parameter
    beta        = ParamCalib(1)             ! discount factor
    G           = ParamCalib(3)             ! government spending
    Debt        = ParamCalib(4)             ! government debt
    ymean       = ParamCalib(5)             ! mean income
    ymeanCalib  = ParamCalib(5)             ! mean income
    
    ! Set tax function parameters
    open(11,  file = 'Input/eqm_results.txt',      status = 'unknown')
    read(11, *) eqm_results
    close(11)
    param_vec_init(1) = eqm_results(3)     ! gamma
    param_vec_init(2) = eqm_results(4)     ! mt
    param_vec_init(3) = eqm_results(5)     ! rt
    omegat = 0d0
    print *, 'Tax parameters initially: gamma, mt, rt = ', param_vec_init
    print *, ''

    ! Guesses for equilibrium objects
    lambda = eqm_results(2)             ! lambda
    r      = eqm_results(1)             ! interest rate
    print*, 'Guesses for r = ', r,' and lambda = ', lambda
    print *, ''
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Compute initial steady state
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! Solve for initial steady state
    print *, 'Computing initial steady-state'
    print *, ''
    
    ! Solve for equilibrium given tax function parameters
    welfare = welfare_eqm(param_vec_init,1)
    
    ! Store eqm objects
    r_init = r                  ! initial equilibrium interest rate
    mt_init = mt                ! initial equilibrium level transfer
    lambda_init = lambda        ! initial equilibrium level tax
    error_vec_init = error_vec  ! initial equilibrium vector of errors
    Kagg_init = Kagg            ! initial equilibrium capital stock

    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Generate grids for grid search over tax systems
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    ! Grids for individual tax-transfer parameters
    mt_grid = LINSPACE_FET(0.10d0, 0.30D0,nmt)           ! transfer level 
    !gamma_grid = LINSPACE_FET(-0.65d0, 0.25D0,ngamma)      ! tax progressivity
    ! Arrange into paramSOB matrix
    paramSOB(1:7,2) = mt_grid(1)
    paramSOB(8:14,2) = mt_grid(2)
    paramSOB(15:21,2) = mt_grid(3)
    paramSOB(22:28,2) = mt_grid(4)
    paramSOB(29:35,2) = mt_grid(5)
    paramSOB(36:42,2) = mt_grid(6)
    paramSOB(43:49,2) = mt_grid(7)
    paramSOB(50:56,2) = mt_grid(8)
    paramSOB(57:63,2) = mt_grid(9)
    paramSOB(64:70,2) = mt_grid(10)
    paramSOB(71:77,2) = mt_grid(11)
    paramSOB(78:84,2) = mt_grid(12)
    paramSOB(85:91,2) = mt_grid(13)
    paramSOB(92:98,2) = mt_grid(14)
    paramSOB(99:105,2) = mt_grid(15)
    paramSOB(106:112,2) = mt_grid(16)
    paramSOB(113:119,2) = mt_grid(17)
    paramSOB(120:126,2) = mt_grid(18)
    paramSOB(127:133,2) = mt_grid(19)
    paramSOB(134:140,2) = mt_grid(20)
    paramSOB(141:147,2) = mt_grid(21)
    
    paramSOB(1:7,1) = LINSPACE_FET(0.24d0, 0.27D0,ngamma)
    paramSOB(8:14,1) = LINSPACE_FET(0.21d0, 0.24D0,ngamma)
    paramSOB(15:21,1) = LINSPACE_FET(0.18d0, 0.21D0,ngamma)
    paramSOB(22:28,1) = LINSPACE_FET(0.155d0, 0.185D0,ngamma)
    paramSOB(29:35,1) = LINSPACE_FET(0.125d0, 0.155D0,ngamma)
    paramSOB(36:42,1) = LINSPACE_FET(0.095d0, 0.125D0,ngamma)
    paramSOB(43:49,1) = LINSPACE_FET(0.07d0, 0.10D0,ngamma)
    paramSOB(50:56,1) = LINSPACE_FET(0.04d0, 0.07D0,ngamma)
    paramSOB(57:63,1) = LINSPACE_FET(0.005d0, 0.035D0,ngamma)
    paramSOB(64:70,1) = LINSPACE_FET(-0.025d0, 0.005D0,ngamma)
    paramSOB(71:77,1) = LINSPACE_FET(-0.06d0, -0.03D0,ngamma)
    paramSOB(78:84,1) = LINSPACE_FET(-0.09d0, -0.06D0,ngamma)
    paramSOB(85:91,1) = LINSPACE_FET(-0.13d0, -0.10D0,ngamma)
    paramSOB(92:98,1) = LINSPACE_FET(-0.165d0, -0.135D0,ngamma)
    paramSOB(99:105,1) = LINSPACE_FET(-0.205d0, -0.175D0,ngamma)
    paramSOB(106:112,1) = LINSPACE_FET(-0.245d0, -0.215D0,ngamma)
    paramSOB(113:119,1) = LINSPACE_FET(-0.285d0, -0.255D0,ngamma)
    paramSOB(120:126,1) = LINSPACE_FET(-0.33d0, -0.30D0,ngamma)
    paramSOB(127:133,1) = LINSPACE_FET(-0.38d0, -0.35D0,ngamma)
    paramSOB(134:140,1) = LINSPACE_FET(-0.43d0, -0.40D0,ngamma)
    paramSOB(141:147,1) = LINSPACE_FET(-0.49d0, -0.46D0,ngamma)
    
    !ix=0
    !do ian = 1, ngamma
    !    do ibn = 1, nmt
    !        ix=ix+1
    !        paramSOB(ix,1) = gamma_grid(ian)
    !        paramSOB(ix,2) = mt_grid(ibn)
    !    enddo
    !enddo
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Compute transitions for points assigned to this program
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    ! How many steady states for each program?
    no_assign = II/no_prog
    
    ! Compute the equilibrium for each point assigned to this program
    loop_start = (seq_no-1)*no_assign+1
    loop_end = (seq_no-1)*no_assign+no_assign
    
    
    do ia = loop_start, loop_end              
       
        ! Show iteration number
        print*, '** iter=', ia
        
        ! Tax-transfer parameter for final steady state
        param_vec_lr(1) = paramSOB(ia,1)    ! tax progressivity
        param_vec_lr(2) = paramSOB(ia,2)    ! transfer level
        param_vec_lr(3) = 0d0               ! phase-out
        
        ! Guesses for eqm objeccts
        r       = 0.02d0
        lambda  = 0.3d0
        
        ! Compute final steady state
        fvalSOB_ss(ia) = welfare_eqm(param_vec_lr,2)
        
        ! Store eqm objects
        r_lr = r                    ! new equilibrium interest rate
        lambda_lr = lambda          ! new equilibrium tax function level
        mt_lr = mt                  ! new equilibrium level transfer
        error_vec_lr  = error_vec   ! new equilibrium vector of errors
     
        ! Show result of this steady state
        print *, ''
        print *, '** done with the steady state with mt = ', mt
        print *, '** Output r, lambda=', r, lambda
        
        ! If steady state converged: transition
        ! 1: asset market; 2: govt BC; 3: VFI; 4: measure
        if (error_vec(1) < 1d-5 .AND. error_vec(2) < 1d-5 .AND. error_vec(3) < 5d-10 .AND. error_vec(4) < 1d-11) then
            
            ! Solving for transition
            print *, ''
            print *, 'Start transition'
            print *, ''
                
            ! Transition
            r_TR = r_lr                                     ! interest rate path
            lambda_TR = lambda_lr                           ! tax level path
            mt_TR = mt_lr                                   ! transfer level path
            wge_TR = (1.0D0-alpha)/(r_TR+delta)             ! wage path
            wge_TR = alpha*(wge_TR**((1d0-alpha)/alpha))    ! wage path
        
            ! Compute Jacobian
            call Compute_JACOB

            ! Compute transition
            call Compute_QNEWTON
        
            
            ! Welfare accounting for transition to new steady state
            fvalSOB(ia) = -sum(mu_TR(:,1)*V_TR(:,1))     
            print *, 'Welfare transition = ', welfare
            print *, ''
        
            ! Store errors
            errors_vec_Sob(ia,:) = error_vec        ! steady state errors
            errors_vec_trans(ia,1,:) = asset_TR     ! transition: deviations from asset market clearing per period
            errors_vec_trans(ia,2,:) = BC_TR        ! transition: deviations from govt BC clearing per period
        
        ! If steady state did not converge: go to next point
        else      
            
            ! Reset eqm objects
            r =  r_init         ! interest rate  
            lambda = 0.3d0      ! tax level
            
            ! Wage implied by interest rate
            wge = (1.0D0-alpha)/(r+delta)
            wge = alpha*(wge**((1D0-alpha)/alpha))
            
            ! Very high value for (negative) welfare
            fvalSOB(ia) = 10000d0
        
            ! Store errors
            errors_vec_Sob(ia,:) = error_vec    ! steady state errors
            errors_vec_trans(ia,1,:) = 1d0      ! transition: deviations from asset market clearing per period
            errors_vec_trans(ia,2,:) = 1d0      ! transition: deviations from govt BC clearing per period
          
        endif
        
        ! Abs. maximum deviations from asset market and govt budget clearing along transition
        error_scalar_trans(ia,1) = maxval(abs(errors_vec_trans(ia,1,:)))  ! asset market 
        error_scalar_trans(ia,2) = maxval(abs(errors_vec_trans(ia,2,:)))  ! BC  
        
        ! Fill info variable
        Info(ia,:) = [ dble(ia), gamma, mt, rt,  error_scalar_trans(ia,1), error_scalar_trans(ia,2), fvalSOB_ss(ia),  fvalSOB(ia),    errors_vec_Sob(ia,3),   errors_vec_Sob(ia,4),  errors_vec_Sob(ia,1),  errors_vec_Sob(ia,2), lambda_lr, r_lr, dble(index_exit), dble(index_it), trans_min_com, mass_trans_min_com, marg_top10_com, marg_avg_com ]
        
        ! Write to file results from info variable
         inquire(file="Output/info.txt", exist=exist)
         if (exist) then
            open(145, file = "Output/info.txt", status="old", position="append", action="write")
         else
            open(145, file="Output/info.txt", status="new", action="write")
         end if
         write(145, '(*(f18.8))') ( Info(ia, isave), isave = 1, N_info)
         close(145)
         enddo
    
end program
